% fig8a_improved.m
% 改进后的吸引域分布图：在 (x0,q0) 平面上，根据改进分类逻辑
% 将吸引子分为：稳定点(1), 周期(2), 混沌(3), 无界(4)

clear; clc; close all;

%% 1. 系统参数设置（根据论文描述调整）
u = 0.1;    % 控制参数
k = 1.8;    % 耦合强度
a = -1;     % memristor参数
b = 15.2;   % memristor参数（此时可能出现周期吸引子，如 period-7）
% 初始条件范围（请根据论文调整）
x0_range = linspace(0, 0.5, 800);    % x0 ∈ [0,0.5]
q0_range = linspace(-0.25, 0.25, 800); % q0 ∈ [-0.25,0.25]

% 预分配结果矩阵，注意这里矩阵的行对应 q0, 列对应 x0
attraction_type = zeros(length(q0_range), length(x0_range));

% 迭代参数
N_total = 2000;     % 总迭代次数
N_discard = 500;    % 舍弃前 200 步
% 注意：阈值的选取需要根据系统实际情况调试
th_stable = 1e-2;   % 稳定点判定标准差阈值
th_period = 1e-3;   % 周期检测阈值
th_unbounded = 1e5; % 无界判定阈值

%% 2. 遍历所有初始条件进行分类
% 注意：我们让矩阵的行对应 q0，列对应 x0
for i = 1:length(x0_range)
    for j = 1:length(q0_range)
        x0_val = x0_range(i);
        q0_val = q0_range(j);
        
        % 迭代系统：调用您的 iterate_nhmm 函数（确保接口正确）
        [x_traj, q_traj] = iterate_nhmm(u, k, a, b, x0_val, q0_val, N_total);
        
        % 取舍弃后轨迹
        x_final = x_traj(N_discard+1:end);
        q_final = q_traj(N_discard+1:end);
        
        % 利用改进分类函数判断吸引子类型
        attraction_type(j,i) = classify(x_final, q_final, th_stable, th_period, th_unbounded);
    end
    fprintf('Progress: %.1f%%\n', i/length(x0_range)*100);
end

%% 3. 绘制吸引域图
figure;
imagesc(x0_range, q0_range, attraction_type);
set(gca,'YDir','normal');  % 保证 q0 从下到上
xlabel('x_0');
ylabel('q_0');
title('吸引域分布');

caxis([1,4]);
colorbar('Ticks',[1,2,3,4], ...
         'TickLabels',{'稳定点','周期','混沌','无界'});

colormap([
    1.0  0.8  0.8;   % (1) 柔和粉色 -> 稳定点
    0.2  0.0  0.4;   % (2) 午夜紫色 -> 周期
    0.6  0.9  1.0;   % (3) 淡青蓝色 -> 混沌
    0.6  0.3  0.1    % (4) 棕色 - 无界
]);
